%
% perc2 - calculation of the fractal dimension of correlated
% bond percolation cluster hulls
%
% Copyright (C) 2009, 2010 Indrek Mandre <indrek(at)mare.ee>
% http://www.mare.ee/indrek/perc2/, http://code.google.com/p/perc2/
% 
% This program is free software; you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation; either version 2 of the License, or
% (at your option) any later version.
% 
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.
% 
% You should have received a copy of the GNU General Public License along
% with this program; if not, write to the Free Software Foundation, Inc.,
% 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
% 

%
% Calculated the sum of the squared residuals as a function
% of the searched exponent. Equivalent to the c++ calc
% application
%
% For octave: following packages required: statistics, optim

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Optionally edit these parameters:

path = 'calc.csv';          % file to read data from
skip = 0 ;                  % how many first data elements to skip
discard = [ ] ;             % which data vectors to discard (1, 2, 3)
exponents = [ ] ;           % add given fixed exponents
range = [ -0.5, 2    ];     % search for data in this range
intervals = 200  ;          % interval points to use
add_index = 0;              % add error to the given column
add_exponent = 2;
add_amplitude = 2;
print = false;

%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

inp = input (['Read data from [' path ']: '], 's');
if numel(inp) ~= 0; path = inp; end
inp = input (['How many elements to skip [' num2str(skip) ']: '], 's');
if numel(inp) ~= 0; skip = str2double(inp); end
inp = input (['Add fixed exponents [' num2str(exponents) ']: '], 's');
if numel(inp) ~= 0
  exponents = [];
  while numel(inp) ~= 0
    [tok, inp] = strtok (inp, '  ,');
    exponents(end+1) = str2double(tok);
  end
end
inp = input (['Discard data columns [' num2str(discard) ']: '], 's');
if numel(inp) ~= 0
  discard = [];
  while numel(inp) ~= 0
    [tok, inp] = strtok (inp, '  ,');
    discard(end+1) = str2double(tok);
  end
end
inp = input (['Search range lower bound [' num2str(range(1)) ']: '], 's');
if numel(inp) ~= 0; range(1) = str2double(inp); end
inp = input (['Search range upper bound [' num2str(range(2)) ']: '], 's');
if numel(inp) ~= 0; range(2) = str2double(inp); end
inp = input (['Number of intervals [' num2str(intervals) ']: '], 's');
if numel(inp) ~= 0; intervals = str2double(inp); end
inp = input (['Add an error to the given column [' num2str(add_index) ']: '], 's');
if numel(inp) ~= 0; add_index = str2double(inp); end
if add_index
  inp = input (['  Error exponent [' num2str(add_exponent) ']: '], 's');
  if numel(inp) ~= 0; add_exponent = str2double(inp); end
  inp = input (['  Error amplitude [' num2str(add_amplitude) ']: '], 's');
  if numel(inp) ~= 0; add_amplitude = str2double(inp); end
end
if print; inp = 'Y/n'; else; inp = 'y/N'; end
inp = input (['Print out data [' inp ']: '], 's');
if numel(inp) ~= 0; print = inp(1) == 'y' || inp(1) == 'Y'; end

fprintf (1, '----------------------------------------------------------\n');

calc = [];
fprintf (1, ['Reading in data from ' path '\n']);
fid = fopen(path, 'rt');
while 1
  tline = fgetl(fid);
  if ~ischar(tline), break, end
  a = [];
  while numel(tline) ~= 0
    [tok, tline] = strtok (tline, ',');
    a(end+1) = str2double(tok);
  end
  if ~isnan(a(1)); calc = cat(1, calc, a); end
end

M = min(calc(:, 2));

if skip ~= 0
  fprintf (1, 'Skipping %d first data elements\n', skip);
  calc = calc(1+skip:end, :);
end

N = size (calc, 1);

if N - M - 1 <= 0
  error ('Not enough data available!');
end

A = calc(:, 1);
X = calc(:, 4:(3+M));
X = X ./ repmat(calc(:, 3), 1, M);
covm = zeros(M, M, N);

data = zeros(M, M);
for i=1:N
  idx = 0;
  for j=1:M
    for k=j:M
      data(j,k) = calc(i,4+M+idx) / calc(i, 3);
      idx = idx + 1;
    end
  end

  for j=1:M
    for k=j:M
      covm(j,k,i) = (data(j,k) - X(i,j)*X(i,k)) / calc(i,3);
      covm(k,j,i) = covm(j,k,i);
    end
  end
end

discard = sort(discard(:),1,'descend');

if numel(discard) ~= 0
  fprintf (1, 'Discarding data column %d\n', discard);
    m = 1:size(X, 2);
  m(discard) = 0;
  X = X(:, find(m));
  covm = covm(find(m),find(m),:);
end

for x=exponents
  fprintf (1, 'Adding a constant exponent %.4f\n', x);
  X = cat(2, X, A .^ x);
  ncovm = zeros (size(covm, 1) + 1, size(covm, 2) + 1, size (covm, 3));
  for i=1:size(covm,3)
    m = cat (2, covm(:,:,i), zeros (size(covm, 1), 1));
    m = cat (1, m, zeros (1, size(m, 2)));
    ncovm(:,:,i) = m;
  end
  covm = ncovm;
end

if add_index
  coeff = add_amplitude * sqrt(covm(add_index, add_index, 1)) / (A(1) ^ add_exponent);
  fprintf (1, 'Adding %.2e*A**%.2f to all rows of column %d (+%.4f for row 1)\n', coeff, add_exponent, add_index, coeff * A(1) ^ add_exponent);
  X(:,add_index) = X(:,add_index) + coeff * A .^ add_exponent;
end

if print
  X
  covm
end

fprintf (1, ['Resulting N=' num2str(size(X,1)) ', M=' num2str(size(X,2)) '\n']);
fprintf (1, ['Searching range [' num2str(range(1)) ', ' num2str(range(2)) '] with ' num2str(intervals) ' intervals..\n']);

fprintf (1, '----------------------------------------------------------\n');

[res, mins] = lsexp (A, X, covm, range, intervals, struct('info', true));

mins

format long e
fprintf (1, 'Writing graph data into ''ocalc.out''\n');
save ('-ascii', 'ocalc.out', 'res');

res;

